Long-time discrete particle effects versus kinetic theory in the self-consistent 

single-wave model 
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The influence of the finite number N of particles cou- 
pled to a monochromatic wave in a collisionless plasma is 
investigated. For growth as well as damping of the wave, 
discrete particle numerical simulations show an TV-dependent 
long time behavior resulting from the dynamics of individ- 
ual particles. This behavior differs from the one due to the 
numerical errors incurred by Vlasov approaches. Trapping 
oscillations are crucial to long time dynamics, as the wave 
oscillations are controlled by the particle distribution inho- 
mogeneities and the pulsating separatrix crossings drive the 
relaxation towards thermal equilibrium. 
PACS numbers: 05.20.Dd (Kinetic theory) 
52.35.Fp (Plasma: electrostatic waves and oscillations) 
52.65.-y (Plasma simulation) 
52.25.Dg (Plasma kinetic equations) 



I. INTRODUCTION 

It is tempting to expect that kinetic equations and 
their numerical simulation provide a fair description of 
the time evolution of systems with long-range or 'global' 
interactions. Yet, physical systems are obviously com- 
posed of a finite (albeit large) number of degrees of free- 
dom. For instance, a plasma is a system of N charged 
particles in interaction. A relevant issue, both for fun- 
damental and computational reasons, is to test whether 
the kinetic description may hide truly physical, finite- N, 
behaviors. 

We address this issue for a typical example where long- 
range interactions come into play, namely wave-particle 
interactions, which are ubiquitous in the physics of hot 
plasmas [0-§- 

Restricting to the electrostatic and collisionless case, 
it has been shown that the interaction of resonant par- 
ticles with Langmuir long-wavelength modes can be de- 
scribed by a Hamiltonian self-consistent model, instead 
of the traditional Vlasov-Poisson kinetic approach. We 
consider the case where particles interact with a single 



mode. In Section II, we present the single-wave Hamilto- 
nian model, whose universal features were recently em- 
phasized in different contexts The essential prop- 
erty of this model is its mean-field nature : resonant 
particles only interact with the wave, which is the only 
effective collective degree of freedom representing longi- 
tudinal oscillations of bulk particles. These bulk par- 
ticles, whose individual characters are absent from the 
physical mechanism at work, namely wave-particle inter- 
action, are then eliminated in our approach. There are 
thus no direct particle-particle interactions in our model. 
There is only a mean-field coupling, which enables one to 
derive rather directly its kinetic limit Q . We use this sin- 
gle wave-particle model to investigate the discrepancies 
between kinetic and finite-V approaches. Numerically, 
the finite-V Hamiltonian microscopic dynamics is com- 
puted through a symplectic code, whereas the numerical 
implementation of the kinetic counterpart of this system 
involves a semi-Lagrangian Vlasov solver. 

In Section III, we focus on long time numerical sim- 
ulations of both kinetic and finite-Y systems for initial 
conditions modeling bump-on-tail beam-plasma instabil- 
ity or damping. We review the well-known results of 
linear theory and show how long time behaviors are in- 
trinsically related to the non-linear regime. Our model 
offers a good alternative to address the recently contro- 
versial issue of the long-time evolution encountered 
in the Vlasov-Poisson system. In Section IV, we analyse 
quantitatively the finite-Y effects. In particular, we show 
how the discrepancies between kinetic and finite-V long 
time behaviors can be related to the presently actively 
investigated topics of the possible inadequacy of Gibbs 
thermodynamics to predict time-asymptotic dynamical 
behaviors in the limit Y — > 00 for long-range systems. 
More precisely, we give hints to the non-commutation of 
limits t — ► 00 and Y — > 00 for the mean-field model. We 
conclude in Section V. 



II. MEAN-FIELD MODEL FOR THE LANGMUIR 
WAVE-PARTICLE INTERACTION 



1 



A. Links between the traditional Vlasov-Poisson 
treatment of Langmuir waves and the self-consistent 
wave-particle Hamiltonian model 

Without collisions, plasma dynamics is dominated by 
collective processes. Langmuir waves and their famil- 
iar Landau damping and growth are a good example 
of these processes, with many applications, e.g. plasma 
heating in fusion devices and laser-plasma interactions. 
For simplicity we focus on the one-dimensional case, rel- 
evant to electrons confined by a strong axial magnetic 
field, and assume that ions act as a neutralizing fixed 
background. The traditional description of the interact- 
ing particles and fields then rests on the (kinetic) coupled 
set of Vlasov-Poisson equations |||| . The current debate 
on the long-time evolution of this system hints that fur- 
ther insight in this fundamental process is still needed 

Our approach to (Langmuir) wave-particle interaction 
complements this usual treatment in that, to capture the 
physical mechanism at work, electrons are partitioned 
in two populations : bulk and tail. The idea behind 
this discrimination is simple : wave-particle interaction 
involves the resonant tail particles whose velocity is close 
to the phase velocity of the wave under consideration. 
These waves are just the collective macroscopic degrees of 
freedom, capturing the longitudinal oscillations of other 
non-resonant, bulk particles, so that these bulk particles 
participate in the effective wave-particle dynamics only 
through the waves. 

Langmuir modes are thus collective oscillations of bulk 
particles, represented by slowly varying complex am- 
plitudes in an envelope approximation. Their interac- 
tion with individual tail particles is described by a self- 
consistent set of Hamiltonian equations [jl 10 111]. These 
already provided an efficient basis for investigating 
the cold beam plasma instability and exploring the non- 
linear regime of the bump-on-tail instability |l3[ ] . Analyt- 
ically, they yield an intuitive and rigorous derivation of 
spontaneous emission and Landau damping of Langmuir 
waves Hj. Besides, as it eliminates the rapid plasma 
oscillation scale u>~ , this self-consistent model offers a 
genuine tool to investigate long-time dynamics. 



B. The single wave model 



Consider an electrostatic potential perturbation 
$(z,t) = 0fe(r) exp i(kz — uj^t) + c.c. (c.c. means com- 
plex conjugate) with complex envelope 4>k, in a one- 
dimensional plasma of length L with periodic boundary 
conditions (and neutralizing background). Wavenum- 
ber k and frequency Uk satisfy a dispersion relation 
e(k,uJk) — 0. The density of N (quasi)-resonant elec- 
trons is cr(z, t) = (nL/N) Yli=i $( z ~ z z( t ))j where n is 
the electron number density and zi is the position at time 
r of electron labeled / (with charge e and mass in). Non- 
resonant electrons contribute only through the definition 
of the dielectric function e, so that 4>k and the z/'s obey 
coupled equations p],|l2[ 



exp[— ikzi + iujki 



~d7 = e k 2 N(de/Ooj k ) ^ 

cPzi iek . 

-—5- = <Pk explikzi - iu k T\ + c.c. 

dr z m 



(1) 



(2) 



with eo the vacuum dielectric constant. We renormalizc 
time and positions to t — ar, xi — kzi — uj^t where a 3 = 
ne 2 /[meo(de/duJk)]- In particular, a — (n/2n p ) 1 / 3 ujp for 
a cold plasma with density n p , plasma frequency uj p , 
and dielectric function e(k,ui) = 1 — uj 2 /oj 2 . With these 
rescaled variables and V — (ek 2 cj)k)/{a 2 m), this system 
defines the self-consistent dynamics (with N + 1 degrees 
of freedom) 



JV 



V = iN- 1 ^2cxp(-ix l ) 



1=1 



Xi = iV exp(zx/) — iV* exp(— ix{) 



(3) 
(4) 



for the coupled evolution (in dimensionless form) of the 
electrons and wave complex amplitude (a star means a 
complex conjugate and ' = d/dt). This evolution derives 
from the Hamiltonian 

JV / 2 \ 
H (x hPl , C, O = E (f - N ~ 1/2 ^ Xl - N-^Ce-^J 

(5) 

where £ = N 1 / 2 V. This system conserves energy TL = H 
and momentum V = YliPl + IC| 2 - An efficient fourth- 
order symplectic integration scheme is used to study this 
Hamiltonian numerically [|l3f . 



We discuss the case of one wave interacting with the 
particles. Though a broad spectrum of unstable waves is 
generally excited when tail particles form a warm beam, 
the single-wave situation can be realized experimentally 
fl5f and allows leaving aside the difficult problem of mode 
coupling mediated by resonant particles Moreover, 
recent studies H|| have stressed the genericness of the 
single-wave model, which we discuss later on. 



C. Kinetic limit and study of finite-iV effects 

As we follow the motion of each particle, we can ad- 
dress the influence of the finite number of particles on the 
long-time behavior of the system. This question is eluded 
by the kinetic Vlasov-Poisson description, and one might 
argue that finite N is analogous to numerical discreti- 
sation in solving kinetic equations. Thus we investigate 



2 



the kinetic limit, N — > oo. As there is no direct particle- 
particle interaction in our model (^), it is possible to 
express in a simple way the N — > oo limit through a 
parallel treatment of the particles and the wave. The 
mean-field coupling between collective (wave) and indi- 
vidual (particles) degrees of freedom enables one to avoid 
the derivation of a full BBGKY hierarchy. 

In the kinetic limit N — > oo, the discrete distribution 
a is thus replaced with a density f(x,p,t), and system 
©"(□) yields the Vlasov-wave system 



dV 
~dt 



—7r=i I exp(—ix)f(x,p,t)dxdp 



(6) 



df df df 
~£it +P dx~ + (iV eXp(ia;) ~ iV * cx PHz)) = 0. (7) 

For initial data approaching a smooth function / as 
N — > oo, the solutions of (^)-(^) have been proved to 
converge to those of the Vlasov-wave system (§)-(0) over 
any finite time interval [Q. This legitimizes our compar- 
ison between finite N and kinetic behaviors. 

The kinetic model @-(0) is integrated numerically by 
a semi-Lagrangian solver, which covers the (x,p) space 
with a rectangular mesh. The function / (interpolated 
by cubic splines) is transported along the characteristic 
lines of the kinetic equation, i.e. along trajectories of the 
original particles ||l7| . Therefore, in addition to the truly 
physical effects of the discrepancies between finite N and 
kinetic systems on long time simulations, we shall also 
compare in this article computational finite grid effects 
of the kinetic solver with the granular aspects of the N- 
particle system Jl8[. 



D. Universal features of the single- wave model 

The single wave model was first formulated as 
a simplified model to treat the instability due to a weak 
cold electron beam in a plasma with fixed ions. For this 
singular case, it was clear that retaining only a single 
Langmuir mode was a good approximation even till some 
primary stage of non-linear saturation. This derivation 
involved natural approximations, but did not a priori 
preserve the Hamiltonian or Lagrangian structure of the 
dynamics (though the latter is recovered in the final equa- 
tions), and a more direct derivation within the Hamil- 
tonian and Lagrangian formalism has been established 

Recently, different studies |2|,|3J have extended the 
regime of application of the single- wave model to a much 
larger class of instabilities and have derived it in a generic 
way in different contexts. 

J.D. Crawford and A. Jayaraman Q studied the colli- 
sionless nonlinear evolution of a weakly unstable mode, 



in the limit of a vanishing growth rate 7 — > + . They de- 
rived in this limit, for a multispecies Vlasov plasma, the 
asymptotic features of the electric field and distribution 
functions. These reveal that the asymptotic electric field 
turns out to be monochromatic (at the wavelength of the 
linear unstable mode) and that the nonresonant parti- 
cles respond to this electric field in an essentially linear 
fashion whereas the resonant particle distribution has a 
much more complicated structure, determined by non- 
linear processes, e.g. particle trapping. That is, starting 
from a much wider class of instabilities than the origi- 
nal single wave model proposed by O'Neil, Winfrey and 
Malmberg |T^| , Crawford and Jayaraman derive asymp- 
totic forms for the electric field and distribution functions 
that precisely feature the assumptions for the single wave 
model. 

D. dcl-Castillo-Negrete j| derived initially the single 
wave picture using matched asymptotic methods to treat 
the resonant and nonresonant particles. In the wider 
context of self-consistent chaotic transport in fluid dy- 
namics, he also showed that the single-wave model pro- 
vided a simplified starting point to study the difficult 
problem of active transport (as opposed to the trans- 
port of passive scalars which do not affect the flow) . Ac- 
tually, the single-wave model captures the essential fea- 
tures of the self-consistent transport of vorticity, i.e. an 
advective field that modifies the flow while being trans- 
ported, through the constraint of a vorticity-velocity cou- 
pling. Self-consistent active transport is a ubiquitous 
phenomenon in geophysical flows or in fusion plasmas 
with the problem of magnetic confinement. 

Finally, we remark that the single wave model has close 
connections, to be further clarified, with systems of cou- 
pled nonlinear oscillators, such as those first studied by 
Kuramoto fLof . The occurrence of a phase transition in 
the regime of Landau damping ^(J , with order parameter 
the mean-field intensity of the wave, is a manifestation 
of these analogies. 

Now we return to the original motivation of this work 
and review wave-beam instability and damping. 



III. WAVE-BEAM INSTABILITY 

A. Linear study 

Let us first study linear instabilities and remark that 
one solution of (||)-(|1) corresponds to vanishing field 
Vq = 0, with particles evenly distributed on a finite set of 
beams with given velocities. Small perturbations of this 
solution have SV = JVbe 7 *, with rate 7 solving Q 



7 = 7r + *7i 



= iN~ 



N 



(l + ipiY 



(8) 
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For a monokinetic beam with velocity U, tffi) reduces to 
7(7 + ill) 2 = i ; the most unstable solution occurs for 
U = 0, with 7 r = \/3/2 and 71 = 1/2. For a warm 
beam with smooth initial distribution /(p) (normalized 
to J fdp =1), the continuous limit of (||) yields 

l = i Jh + ip)- 2 f(p)d P . (9) 

For a sufficiently broad distribution (|/'(0)| <C 1), we 
obtain |7 r |7 r = 7r7r/'(— 7i), where /' = df /dp, and, for 
|/"(0)| < 7T- 1 , one finds 71 w 7T7 r /"(0). Except for the 
trivial solution given by 7 r = 0, one easily checks that 
other solutions can only exist for a positive slope /'(0). 
Then the perturbation is unstable as the time behavior 
of SV is controlled by the eigenvalue 7 with positive real 
part, i.e. with growth rate 7 r w 7l = 7r/'(0) > 0. For 
negative slope, one recovers the linear Landau damping 
paradox : the observed decay rate 7l = tt/'(0) < 
is not associated with genuine eigenvalues, but with 
phase mixing of eigenmodes @,^-|||. This is a direct 
consequence of the Hamiltonian nature of the dynamics 




p < 3.96 and f(p) — otherwise, with a = 11.89 and 
pi = 15.85, 

(Hi) a truncated Lorentzian (TL) with positive slope 
f(p) = (c/tt)/[(p - pi) 2 + a 2 ] if -7.42 < p < 3.18 and 
/(p) = otherwise, with a = 2.12 and p\ — 1.06. 

For the growing case, nonlinearities result from the 
growth of the wave intensity. For the damping case, the 
linear description introduces time secularities which ulti- 
mately may cause the linear theory to break down. The 
ultimate evolution is intrinsically nonlinear, not only if 
the initial field amplitude is large, as in O'Neil's seminal 
trapping picture for one wave j8j, but also if one con- 
siders the system evolution over time scales of the order 
of the trapping time (which may be large if the initial 
wave amplitude is small). The question of the long-time 
fate of plasma wave amplitude is thus far from trivial 
H . Though some simulations Q suggest that nonlinear 
plasma waves eventually approach a Bernstein-Greene- 
Kruskal steady state (2fJ instead of a Landau vanishing 
field, the answer should rather strongly depend on ini- 
tial conditions ||. Our iV-particle, 1-wave system is the 
simplest model to test these ideas. 



B. Nonlinear regime 

This linear analysis generally fails to give the large 
time behavior. This is obvious for the unstable case as 
non-linear effects are no longer negligible when the wave 
intensity grows so that the bounce frequency u)b(t) — 
y/2\V(i)\ of trapped particles in the wave becomes of the 
order of the linear growth rate j T . 

We used the monokinetic case as a testbed ^l],^3| . As 
seen in Fig. [j], finite-iV simulations show that the un- 
stable solution grows as predicted, until it saturates to 
a limit-cycle-like behavior where the trapping frequency 
u>b(t) oscillates between 1.27 r and 27,.. In this regime, 
some of the initially monokinetic particles have been 
scattered rather uniformly over the chaotic domain, in 
and around the pulsating resonance, while others form a 
trapped bunch inside this resonance (away from the sep- 
aratrix) as observed in Fig. ^| p3|] . This dynamics is quite 
well described by effective Hamiltonians with few degrees 
of freedom (n],^l| . Note that it cannot be easily followed 
by a numerical Vlasov solver, as the initial beam has a 
singular velocity distribution function. 

In this article, we discuss the large time behavior of 
the warm beam case, with f'(po) ^ at the wave nomi- 
nal velocity po = 0. Figure ^ displays three distribution 
functions (in dimcnsionlcss form) with similar velocity 
width (here c normalizes J fdp — 1 in each case) : 

(i) a function (CD) giving the same decay rate for all 
phase velocities, /(p) = c~ac/(pi—p) if —3.96 < p < 3.96 
and /(p) = otherwise, with a = 11.89 and p\ — 15.85, 

(ii) a function (CG) giving a constant growth rate for all 
phase velocities [Q, /(p) = c— ac/{p\ + p) if —3.96 < 



IV. FINITE TV-EFFECTS AND KINETIC 
TREATMENT IN LONG-TIME DYNAMICS 

First of all let us mention that for finite N, the particles 
are initially distributed in (x,p) so that their distribution 
approaches / smoothly in the large N limit ]l|, p^| , pT[ . 

A. The damping case 

A thermodynamical analysis |EQ] predicts that, for a 
warm beam (i.e. if the velocity distribution / has a large 
width, as in Fig. [jj) and small enough initial wave am- 
plitude, Wb asymptotically scales as TV -1 / 2 in the limit 
N — ► 00. Figure || shows the evolution of a small ampli- 
tude wave launched in the beam. The iV-particle system 
(curve N) and the kinetic system (curve V) initially damp 
exponentially as predicted by perturbation theory p4[ , 
for a time of the order of |7l| _1 - After that phase-mixing 
time, trapping induces nonlinear evolution and both sys- 
tems evolve differently. For the N -particle system, the 
wave grows to a thermal level that scales as TV" 1 / 2 , corre- 
sponding to a balance between damping and spontaneous 
emission |l~4|,|2(|. For the kinetic system, initial Landau 
damping is followed by slowly damped trapping oscilla- 
tions around a mean value ; this mean value also decays 
to zero, at a rate which decreases for refined mesh size. 
Figure [| thus reveals that finite N and Vlasov behav- 
iors can considerably diverge as spontaneous emission is 
taken into account. 

Figure || represents the time evolution of the wave am- 
plitude for different values of N. It clearly shows how 
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finite-TV wave evolutions depart progressively from the 
N — ► oo curve (the later for larger N). One should also 
notice, from Figs [| and [3] and for sufficiently large N, 
the onset of nonlinear effects at large time. In spite of 
the smallness of the initial values of the wave amplitude, 
nonlinear effects (through trapping) eventually come into 
play and stop Landau exponential decay, marking the 
beginning of a different dynamical regime for which the 
decay is far slower. 



B. The single wave warm beam instability 

Now consider a warm beam with a velocity distribu- 
tion with a positive slope at po — 0. Line Nl (resp. 
N2) of Fig. H displays m(cJb(i)/7r) versus time in a nu- 
merical integration of @-(f|) for a CG distribution with 
N = 128000 (resp. 512000) and -y T = 0.08. Line VI (resp. 
V2) shows the evolution of ln(wb(i)/7r) versus 7 r i in nu- 
merical integration of the kinetic system for the same 
initial conditions with a 32 x 128 (resp. 256 x 1024) grid 
in (x,p) space. All four curves exhibit the same initial 
exponential growth of linear theory with less than 1% er- 
ror on the growth rate. Saturation occurs for u>b/7r ~ 3.1 
P). Lines Nl and VI do not superpose beyond the first 
trapping oscillation after saturation. Note that, in our 
system, oscillating saturation cannot be related to ex- 
citation of sideband or satellite Langmuir waves as our 
single- wave Hamiltonian does not allow for any spectrum 
of waves. 

Beyond the first trapping oscillation, kinetic simula- 
tions exhibit a second growth at a rate controlled by mesh 
size. Line V2 suggests that a kinetic approach would pre- 
dict a level close to the trapping saturation level on a time 
scale awarded by reasonable integration time. This level 
is fairly below the thermodynamic level Vth predicted by 
a Gibbsian approach |2C|] . Such pathological relaxation 
properties in the N — > oo limit seem common to mean- 
field long-range models ||(|. Both kinetic simulations 
also exhibit a strong damping of trapping oscillations, 
which disappear after a few oscillations, whereas finitc- 
N simulations show persistent trapping oscillations. 

One could expect that finite-TV effects would mainly 
damp these oscillations, so that the wave amplitude 
reaches a plateau. However their amplitude does not 
depend on the number of particles N, which shows that 
they are not an artifact due to 'poor accuracy' of finite- TV 
simulations. Moreover the wave amplitude slowly grows 
furl her. wlinvus the velocity distribution function flat- 



tens over wider intervals of velocity |18 



This second growth after the first trapping saturation 
depends on the shape of the initial distribution function. 
In Fig. ||(b), curve N2 is the same as in Fig. ^(a) but 
computed over a longer duration, and curve N3 corre- 
sponds to N = 64000 with the TL distribution of Fig. |. 



Although curve N3 corresponds to 8 times fewer parti- 
cles than curve N2, the final level reached at the end of 
the simulation is lower. In the second growth regime, 
particles are transported further in velocity, so that the 
plateau in f(p) broadens with time. As will be clearly 
shown in Sec. D, this spreading is due to separatrix cross- 
ing, i.e. successive trapping and detrapping by the wave 
p3[ . As the resonance width of the wave separatrix 
grows, the wave can trap particles with initial velocity 
further away from its phase velocity. Noting that the TL 
distribution decays while the CG distribution still grows 
for v > 0.05, we see that, with TL, fewer particles can 
give momentum to the wave when being trapped (as V 
is conserved) ; hence the second growth is slower for the 
TL distribution. 

We followed the evolution of the wave amplitude of 
curve N3 of Fig. ||(b) up to j T t — 1750. Starting from the 
first trapping saturation level, equal to 40% of the ther- 
modynamic level Vth, we observe persistent amplitude 
fluctuations with a growth rate that slowly decreases as 
we reach 0.78Vth at the end of the computation. 

Line N4 of Fig. |^(b) corresponds to the TL distribu- 
tion with 2048000 particles and shows persistent oscil- 
lations with approximately the same amplitude as for 
N = 64000. 



C. Trapping oscillations 

Let us show that the occurrence of trapping oscillations 
with non-vanishing amplitude follows from the existence 
of spatial inhomogeneities. For this purpose, introduce 
the complex field 



1 N 

exp(ixi) = exp(ia 



i=i 



(10) 



In the kinetic Vlasov limit, it is the first Fourier compo- 
nent of the spatial density 



M (oo) = / exp(ix)f(p,x,t)dpdx. 



(11) 



A spatially homogeneous phase space with independent 
particles corresponds obviously to a |AfW| scaling as 
N^ 1 / 2 i.e. to a vanishing |M(°°)|. From (||), and drop- 
ping the superscript (N), it follows that 



V = iM*. 

Putting 

2V = 2 \V\ exp(-i6) = u\ exp( 
one obtains from ( [l2| ) 
d\V\ 



dt 



|M|sin(a - 6). 



(12) 
(13) 

(14) 
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Moreover, if the wave amplitude displays a sinusoidal 
temporal evolution (e.g. due to trapping oscillations, 
then toi = u>ho) such that 

\V\(t) = ^l + AVsm(uj 1 t) (15) 

with Ay > the amplitude of oscillations and Wbo the 
(quadratic) average bounce pulsation, one obtains 

d\V\ 



dt 



u>i AV cos(ujit) 



(16) 



so that taking the time average, denoted by (-) t , of the 
square of both (Ful) and (Il6|) one gets 



(\M?Y t 



1/2 



(17) 



provided that a — 9 has a uniform distribution (e.g. if 
a — 9 LOit for some lo-i). This simplified model, sup- 
posing only a harmonic oscillation for |V|, shows that 
the amplitude of the wave oscillations depends directly 
on the occurrence of (x,p)-space inhomogeneities. Ac- 
tually, for a homogeneous phase space, ( |l~7| ) implies that 
Ay scales at most as TV^ 1 / 2 and vanishes in the kinetic 
limit. 

However, consider now the case where a — 9 is the 
barycentric position, in the reference frame of the wave, 
of an inhomogeneity (clump) composed of a finite fraction 
\M\ of the particles. If this clump is sufficiently close to 
the bottom of the potential well, then a(t) — 9(t) is small 
and (|l7j ) can be reduced to 



Ay = uj- 1 \M\(2(a- 



(18) 



? 1/2 

where ((a — #) 2 ) t is the mean distance of the clump 
from the center of the resonance cat's eye. Note that 
Ay may be small even if M is large, provided that the 
clump stays close to the bottom of the well [^Tj . 

A striking example is given by the cold beam-wave in- 
stability jll| , for which a macroscopic fraction of the par- 
ticles belong to a so-called macroparticle that oscillates 
near the bottom (elliptic point) of the potential well and 
drives the wave amplitude oscillations as shown in Fig. ||. 
Another illustration of (|l^) is given by Figs 2 and 4 of 
ref. H . In Vlasov simulations, as the spatial resolution 
N x is increased, one observes a more refined phase space 
that reveals some heterogeneities for the larger value of 
N x . This simultaneously goes with pronounced trapping 
oscillations at large N x while the rough resolution, which 
has smeared out the thin filamcntation in the vortices, is 
associated to a flat, constant amplitude in time. 



D. Finite-iV effects, trapping oscillations and 
relaxation towards thermal equilibrium through 
chaos 

We now discuss the actual process by which the wave- 
particle system relaxes. For this purpose, one can get a 



flavor of the stochasticity (strong or weak depending on 
resonance overlapping or not) that a test particle would 
encounter in (x,p) space at different stages of the evolu- 
tion. 

In the equation of motion (|^) of any particle in the 
system (||) 



x = -wgW sin ( x ~ #(*)) 



(19) 



the time dependence of the wave modulates the force on 
the test particle. If V is approximately periodic over 
a time interval of length T, and if its period is large 
compared with the trapped particle bouncing period, the 
pulsations of the separatrix generate strong chaos in the 
domain of (x,p) space it sweeps ^8|. In the present case, 
the period of V is basically the trapped particles' bounc- 
ing period, which makes the 'slow chaos' approximation 
rather crude. 

The chaos generated by a pulsating separatrix can 
also be characterized using the Fourier decomposition 
of the wave V(t) = J2™ n =-ca |vn|exp(iw n i + Xn), with 
uj n = 2nn/T. Then the test particle experiences a force 
deriving from an effective many-wave field (though the 
Hamiltonian @ involves a single wave), 



-2 ^ \ v n\ sin(a; + w n t - Xn) 



(20) 



and the overlaps between resonances in this force field 
cause the particle to move chaotically |2!| . 

We computed the Fourier decomposition of u>b(t) for 
the warm beam with initial distribution TL over two 
different time intervals, one (window T\) just after the 
nonlinear saturation with 38 < 71^ < 65 and the other 
one (window T2) far later with 250 < 71^ < 302, for 
N = 48000 and N = 768000 particles. The low fre- 
quency bias induced by the slow growth of the amplitude 
evolution has been removed by subtracting the quadratic 
best fit of the wave evolution. The Fourier decomposition 
(with 2tt /T <^ ujbo and T = T\ or T2) reads 



u)?(t) = a ( 1 + — cos(w„t + ip n ) 



where 



a 



J b0 



T 



to+T 



ul(t)dt 



(21) 



(22) 



with 7l£o = 38 or 250. Figs ^ and || show coefficients 
o n /oo as a function of the frequency normalized by WbOi 
the fundamental n = being removed. Such a nor- 
malized form enables a direct comparison between the 
figures (irrespective of N and T). Below we use co- 
efficients cihn = On', where n' is the index such that 
lu h i — n'2TT /T — nuj^Q. 

During the first stage T\ (Fig. [7]), the behavior appears 
almost entirely driven for large N by a narrow spectrum 



6 



of frequencies of the order of the average trapping fre- 
quency, although, for N = 48000, other Fourier compo- 
nents are excited. 

Two types of chaos must be distinguished. The sim- 
plest one is related to the classical overlap between the 
resonances of two waves p9| propagating at velocities 
TOCJbo and ra^bo- The corresponding chaos is 'fast', 
namely the frequencies of the relevant resonances are 
larger than cjbo- The stochasticity parameter Q esti- 
mated as 

V^flbri + V^ a bm _ y/^bn/ a + \/2abm/ ^0 
s m,n — j j — 7~ j 

\m — n\ujhQ \tn — n\ 

(23) 

is small for all m 7^ 0, n 7^ 0. This means that the 
effective many-wave field does not generate strong chaos 
in velocity ranges away from the wave nominal velocity. 

Similarly, one may search for chaos induced by the 
overlap between the component in ( po|) with phase ve- 
locity Wbn > 2wbo, and the central resonance, with 
phase velocity 0. Their resonance overlap |3(| parameter 
so, ri = (2-y/eto + V2abn)/|?i^bo| ~ 2/n is not large enough 
in our system to induce large scale chaos, transporting 
particles from the neighborhood of the wave nominal ve- 
locity (with its natural width 2wbo) to the neighborhood 
of waves with significantly different phase velocities. 

The second type of chaos is related to the resonant 
forcing by Fourier components with phase velocity Wbn ^ 
2wb0i on t ne periodic motion induced by the main field 
component. This process applies in the vicinity of the 
separatrix of the main resonance and must be discussed 
using action-angle variables. As recalled in the appendix, 
the corresponding pulsations are smaller than 2wbo for 
untrapped particles, and smaller than Wbo for trapped 
particles (and for untrapped ones very close to the sepa- 
ratrix). Particles experiencing this chaos easily cross the 
pulsating separatrix, i.e. change between trapped and 
untrapped motion. 

The Fourier spectra in Fig. show that the second 
process is active, since only to < Wbo lead to significant 
amplitudes a^ n . This supports the analysis of chaos in 
our system as 'slow chaos' due to the pulsating resonance 
p8|| . However, corresponding amplitudes a^n are much 
smaller than clq. Therefore, the pulsating resonance does 
not sweep the vicinity of the bottom of the wave potential 
well, which allows the particles close to the elliptic point 
to move quite regularly, forming a macroparticle like in 
Fig. ||. This bottom of the well may be separated from 
the surrounding chaotic domain (swept by the separatrix) 
by KAM surfaces in (x,p,t) space Q. 

The more chaotic behaviour is expected in the case 
where the additional peaks, close to uJbOi 2c^bo and 3o>b0j 
are more important. As this is the case N — 48000, our 
analysis is consistent with the faster transport in velocity 



observed for the smaller values of N and thus with the 
more rapid thermalization in this case. 

Considering the later time interval T 2 (Fig. ||), it is 
striking to note that the differences between both spec- 
tra are strongly reduced, so that one can estimate that 
a test particle feels a similar amount of stochasticity for 
both values of N. Moreover, the relative intensity of the 
secondary resonances is smaller than in the time inter- 
val T\ , so that the separatrix pulsations will be relatively 
smaller. As a result, the chaotic growth of the wave also 
slows down. By contrast, the Vlasov simulations, with 
their induction of coarse graining, even prevent the evo- 
lution of the system towards this second stage. Actually, 
these codes are unable to capture the spatial details of 
the phase space filamentation which are under the scale 
of the mesh size, so that, after a certain stage, they arti- 
ficially make the system almost integrable. 

E. (i,p)-space structures, smoothing and numerical 
entropy production in the Vlasov model 

Finally, our numerical observations along with equa- 
tion ( |l7|) enable one to describe properly the finite-TV 
effects on the thermalization process. Due to the ini- 
tial asymmetry in velocity space (the initial distribution 
function has a positive slope around the wave phase ve- 
locity), heterogeneities will exist in the wave resonance 
during the first nonlinear oscillations. If N is increased, 
the (weak) chaos available to thermalize the system dis- 
appears, and the system is driven by the filamentary 
structure which develops inside the wave potential well. 
The bottom of the well is an elliptic point for the model, 
around which the correlations between trajectories may 
decay algebraically |3l]] , so that the filaments get thinner 
and more entangled as time goes on. 

The filamentation described here mixes particles in the 
neighborhood of the wave resonance. Particles with ve- 
locities far away from the resonance are weakly sensi- 
tive to the resonance, which they essentially average off. 
Therefore, as N becomes very large, the system long-time 
evolution is dependent on initial conditions in the neigh- 
borhood of the resonance, and thermodynamical conclu- 
sions (relying on ergodicity in the energy surface and ba- 
sic conservation laws) no longer apply. 

Now, to what extent do actual Vlasov simulations 
reproduce either the kinetic equation evolution or the 
finite-iV evolution ? Crudely speaking, kinetic simula- 
tions start to induce a bias with respect to the Vlasov 
equation, because of phase space averaging, when the 
distribution function exhibits structures on scales below 
their grid resolution (and this is bound to happen). One 
might expect that finer grids would enable one to de- 
scribe more precisely the long-time evolution of the sys- 
tem. However, refined grids also reduce the coarseness of 
the particle distribution in (x,p) space, and inhibit the 
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noisy separatrix crossings. This in turn inhibits the wave 
second growth, which results from the coarseness of the 
particle distribution for an actual finite value of N. 

Our observations indicate that kinetic models are too 
idealized, in comparison with finite N, and do not con- 
tain all the intricate behavior displayed by a discrete par- 
ticles system. This can also be tested within kinetic the- 
ory. In particular, whereas the kinetic equation analyti- 
cally preserves integral functionals of / like the 2-entropy 
<?2(t) = /(I — f)fdxdp, numerical schemes increase S2 
significantly, as shown in Fig. |9|, when constant-/ con- 
tours form filaments in (x,p)-space on scales below the 
grid mesh. This filamentation (due to anisochronism of 
nonlinear trapped motion, or shear in (x,p)-space) is 
smoothed by numerical partial differential equation in- 
tegrators, while TV-body dynamics follows the particles 
more realistically, sustaining the trapping oscillations. 

This inability of the semi-Lagrangian scheme in re- 
producing the filamentation over small scales is shared 
by other Vlasov schemes. In particular, particle-in-cell 
(PIC) schemes explicitly replace the actual N particles 
of the initial dynamics by effective particles, which are 
redefined smoothly at each time step. PIC schemes 
have a further disadvantage in comparison with semi- 
Lagrangian schemes. They put the numerical effort on 
cells according to their populations ; by contrast, semi- 
Lagrangian schemes ensure a similar accuracy for the 
poorly populated (a;,p)-space domains as for the highly 
populated ones jl7|], so that they describe frontiers in / 
more sharply. 

This discussion shows that the 'irreversible' growth of 
the wave is not related to the entropy production of this 
sub-grid-scale filamentation but to the chaotic trapping- 
detrapping process (of which some small-scale structures 
are by-products). Although the smoothing may appear 
as a minor nuisance in the chaotic (x,p) regions, it does 
actually force a distinctly different evolution in the long 
term, and refining the mesh does not prevent this. 

V. COMMENTS AND CONCLUSION 

In summary, dealing with the basic propagation of 
a single electrostatic wave in a warm plasma, we pre- 
sented finite-iV effects which do not merely result from 
numerical errors and are missed in a kinetic simulation 
approach. Their understanding depends crucially on de- 
scribing the particle dynamics in (x,p) space. The sen- 
sitive dependence of the microscopic chaotic evolution 
to the fine structure of the initial particle distribution 
in (x,p) space pi) implies that the limits t — > 00 and 
N —* 00 do not commute. The driving process in the 
system evolution is separatrix crossing, which requires 
a geometric approach to the system dynamics. Further 
work in this direction will also shed new light on the 
foundations of frequently used approximations, such as 



replacing original dynamics by coupled stochastic 

equations, in which particles undergo noisy transport. 
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A. Appendix : Periodic orbits of the pendulum and 
resonance overlap 

The pendulum equations reduce to the normal form 
( |l9| ) with fixed parameters 9 = and ujho- The orbits of 
the pendulum in (x,p) space (with x modulo 2tt) are : 
fixed points (0, 0) and (tt, 0), the two branches of the sep- 
aratrix and the three types of periodic orbits. They are 
parametrized by the energy E = p 2 /2+w^ (l — cos x) > 
and (if untrapped) by the sign of p. At the unstable fixed 
point and on the separatrix, E = 2cl> 2 . 

The separatrix branches are limits of periodic or- 
bits with periods going to 00. Their equation p± = 
±2ajbo cos(x/2) implies that the velocity ranges in 

[— 2^bo, 2wboJ- 

Circulating orbits have periods decreasing for increas- 
ing energy E = 2uj 2 k~ 2 , where < k < 1. Their period 
is T(k) = 2fcw b ~ 1 K(fc 2 ), with the complete elliptic inte- 
gral K, and k — > I on approaching the separatrix. Period 
T = 27njJ b0 1 occurs for k s» 0.99, i.e. extremely close 
to the separatrix, with E m 2.04 w 2 . The next strong 
resonance, with T = hlj^q , occurs for k rj 0.8, i.e. for 
particles with E w 3.17 a^o and 1.5 c^bo — \p\ — 2.5 Ubo- 

Trapped orbits have periods T(k) = 4w b0 1 K(fc 2 ), with 
< k < 1. The period is larger than 2n/ojbo and in- 
creases with the energy E = 2ui 2 iQ k 2 . The velocity is in 
the range [-y/2E,y/2E]. 

To apply the resonance overlap picture, one considers 
only the relative velocity of the two waves whose reso- 
nances overlap. As the range associated with the main 
cat's eye is [— 2ojbo, 2ojbo], the classical resonance overlap 
picture makes sense only for waves with relative velocity 
larger than 2ajbo with respect to the principal one. 
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FIG. 4. Time evolution of ln(wb(£)/|7L,|) for a CD velocity 
distribution and initial wave amplitude below thermal level : 
(N) JV-particles system with N = 32000, (V) kinetic scheme 
with 32 x 512 (x,p) grid. Inset : short-time evolution. 
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FIG. 5. Time evolution of ln(wb(t)/|7L,|) for an initial CD 
velocity distribution and different values of N. When neces- 
sary for readability, curves were smoothed through a sliding 
window average. 
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FIG. 6. Time evolution of In (ojb(i)/7r)- (a) CG initial dis- 
tribution : kinetic scheme with (VI) 32 x 128, (V2) 256 x 1024 
(x,p) grid ; iV-particles system with (Nl) N = 128000, 
(N2) N = 512000 ; (b) Comparison of CG (N2) with TL 
initial distribution for (N3) N = 64000, (N4) N = 2048000. 




FIG. 8. Coefficients a n /ao (n > 1) of the Fourier decom- 
position of amplitude u!^(t) during the second time window 
250 < 7l£ < 302 as a function of normalized frequency, for 
(a) N = 48000 and for (b) N = 768000 particles. 
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FIG. 7. Coefficients a n /ao (n > 1) of the Fourier de- 
composition of amplitude u£(t) during the first time window 
38 < 7 L t < 65 for N = 48000 (solid curve) and N = 768000 
(dashed curve) particles, as a function of normalized fre- 
quency. 
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FIG. 9. Evolution of the 2-entropy S 2 = J (I - f)fdxdp 
as a function of time 7 r t in the Vlasov simulations for several 
N x x grids. The departure from zero indicates that the 
damping of trapping oscillations is spurious. 
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